home *** CD-ROM | disk | FTP | other *** search
/ Amiga Tools 3 / Amiga Tools 3.iso / grafik / raytracing / rayshade-4.0.6.3 / libray / libobj / mountain.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-09  |  26.1 KB  |  840 lines
  1. /*
  2.  * mountain.c
  3.  *
  4.  * Peter Janssen
  5.  *
  6.  */
  7.  
  8. #include "geom.h"
  9. #include "mountain.h"
  10. #include "allocmatrix.h"
  11. #include <math.h>
  12. #include <stdio.h>
  13. #include <stdlib.h>
  14.  
  15.  
  16. static void MountainEvolve();
  17. static void MountainSetupLevelList();
  18. static void MountainSetupCache();
  19. static void MountainSetupCalcList();
  20. int MountainBoundsIntersect();
  21. int CheckGrid();
  22. int IntersectTriangle();
  23.  
  24. #define    x2grid(m, x) (((x) - m->D2Bounds[0][0]) * m->InvXSize)
  25. #define    y2grid(m, y) (((y) - m->D2Bounds[0][1]) * m->InvYSize)
  26.  
  27. #define    grid2x(m, x) ((x) * m->XSize + m->D2Bounds[0][0])
  28. #define    grid2y(m, y) ((y) * m->YSize + m->D2Bounds[0][1])
  29.  
  30. #define OutOf2DBounds(p, b)                                 \
  31.                ((p).x < (b)[0][0] || (p).x > (b)[1][0] ||   \
  32.                 (p).y < (b)[0][1] || (p).y > (b)[1][1] )
  33.  
  34. static Methods *iMountainMethods = NULL; 
  35. static char mountainName[] = "mountain";
  36. unsigned long MountainTests, MountainHits;
  37. static int WhichCache[3][3] = { {TRUE, TRUE, FALSE},
  38.                                 {FALSE, FALSE, FALSE},
  39.                                 {FALSE, TRUE, TRUE}};
  40.  
  41.  
  42. int MyPower(a, b)
  43. int   a, b;
  44. {
  45. int   result, i;
  46.  
  47.    for (result = 1, i = 1; i <= b; i++) result *= a;
  48.    return result;
  49. }
  50.  
  51. /* Random Generator */
  52.  
  53. #define rand_m  (unsigned long)2147483647
  54. #define rand_q  (unsigned long)127773
  55.  
  56. #define rand_a (unsigned int)16807
  57. #define rand_r (unsigned int)2836
  58.  
  59. /*
  60.  * F(z) = (az)%m
  61.  *      = az-m(az/m)
  62.  *
  63.  * F(z)  = G(z)+mT(z)
  64.  * G(z)  = a(z%q)- r(z/q)
  65.  * T(z)  = (z/q) - (az/m)
  66.  *
  67.  * F(z)  = a(z%q)- rz/q+ m((z/q) - a(z/m))
  68.  *       = a(z%q)- rz/q+ m(z/q) - az
  69. */
  70.  
  71. static long rand_seed;
  72.  
  73. static long rand_number()
  74. {
  75. long    lo, hi, test;
  76.  
  77.     hi   = rand_seed/rand_q;
  78.     lo   = rand_seed%rand_q;
  79.  
  80.     test = rand_a*lo - rand_r*hi;
  81.  
  82.     if (test > 0)
  83.         rand_seed = test;
  84.     else
  85.         rand_seed = test + rand_m;
  86.  
  87.     return rand_seed;
  88. }
  89.  
  90. double GAUSS(s, x, y, level)
  91. int  s, x, y, level;
  92. {
  93.    long    x1, x2;
  94.    double  r, v1, v2;
  95.    static double inv_rand_m;
  96.    inv_rand_m = 1.0 / (double) rand_m;
  97.  
  98.    rand_seed = s + (y<<15) + x + (level<<8);
  99.    do {
  100.        x1 = rand_number();
  101.        x2 = rand_number();
  102.        v1 = (((double)x1 * inv_rand_m) * 2) - 1;
  103.        v2 = (((double)x2 * inv_rand_m) * 2) - 1;
  104.        r = v1*v1 + v2*v2;
  105.    } while (r >= 1.0);
  106.    r = sqrt(-2.0 * log(r) / r);
  107.    return v1*r;
  108. }
  109.  
  110.  
  111. /* 
  112.  * Create MountainObject and return Reference to it
  113.  */
  114. Mountain *MountainCreate(p1, p2, number, dimension, evolvelevel,
  115.                          iterationlevel, scale, startgrid)
  116. Vec2d            p1, p2;
  117. int              number, evolvelevel, iterationlevel;
  118. Float            dimension, scale;
  119. StartGridStruct *startgrid;
  120. {
  121.         Mountain *mountain;
  122.         Vec2d     size;
  123.         int       gridmax;
  124.  
  125.         mountain = (Mountain *)share_malloc(sizeof(Mountain));
  126.         size.u = fabs(p1.u - p2.u);
  127.         size.v = fabs(p1.v - p2.v);
  128.         if (equal(size.u, 0.) || equal(size.v, 0.)) {
  129.             RLerror(RL_ADVISE, "Degenerated startgrid.\n");
  130.             return (Mountain *) NULL;
  131.         }
  132.         mountain->D2Bounds[LOW][X] = min(p1.u, p2.u);
  133.         mountain->D2Bounds[HIGH][X] = max(p1.u, p2.u);
  134.         mountain->D2Bounds[LOW][Y] = min(p1.v, p2.v);
  135.         mountain->D2Bounds[HIGH][Y] = max(p1.v, p2.v);
  136.  
  137.         if (number < 0 || number > 65565) {
  138.             RLerror(RL_WARN, "Illegal MountainNumber : %d", number);
  139.             number = 1000;
  140.         }
  141.         if (dimension < 0.5 || dimension > 3.5) {
  142.             RLerror(RL_WARN, "Illegal FractalDimension : %f", dimension);
  143.             dimension = 2.5;
  144.         }
  145.         if (iterationlevel < 0) {
  146.             RLerror(RL_WARN, "Illegal IterationLevel : %d", iterationlevel);
  147.             iterationlevel = 2;
  148.         }
  149.         if (evolvelevel < 0 || evolvelevel > iterationlevel) {
  150.             RLerror(RL_WARN, "Illegal EvolveLevel : %d", evolvelevel);
  151.             evolvelevel = 0;
  152.         }
  153.         if (scale < 0) {
  154.             RLerror(RL_WARN, "Illegal ScaleFactor : %f", scale);
  155.             scale = 1.0;
  156.         }
  157.  
  158.         gridmax  = MyPower(2, iterationlevel) * (startgrid->Size - 1) + 1;
  159.         mountain->XSize    = size.u / gridmax;
  160.         mountain->YSize    = size.v / gridmax;
  161.         mountain->InvXSize = 1.0 / mountain->XSize;
  162.         mountain->InvYSize = 1.0 / mountain->YSize;
  163.  
  164.         mountain->Number           = number;
  165.         mountain->IterationLevel   = iterationlevel;
  166.         mountain->EvolveLevel      = evolvelevel;
  167.         mountain->FractalDimension = dimension;
  168.         mountain->Scale            = scale;
  169.         mountain->StartGridSize    = startgrid->Size;
  170.         mountain->StartGridData    = startgrid->Data;
  171.         mountain->GridMax          = gridmax;
  172.  
  173.  
  174.         MountainEvolve(mountain);
  175.         MountainSetupLevelList(mountain);
  176.         MountainSetupCache(mountain);
  177.         MountainSetupCalcList(mountain);
  178.  
  179.         return mountain;
  180. }
  181.  
  182. /*
  183.  * Calculating MountainData up till iterationlevel 'EvolveLevel'
  184.  */
  185. static void MountainEvolve(m)
  186. Mountain *m;
  187. {
  188.    Float      Scale;
  189.    Float    **OldMountain, **NewMountain;
  190.    int        OldSize, NewSize;
  191.    int        XFar, XShort, YFar, YShort;
  192.    int        i, x, y;
  193.  
  194.    Scale         = m->Scale;
  195.    OldMountain   = m->StartGridData;
  196.    OldSize       = m->StartGridSize;
  197.  
  198.    for (i = 1; i <= m->EvolveLevel; i++) {
  199.        NewSize = 2 * OldSize - 1;
  200.        Scale = Scale * pow((OldSize / (Float) NewSize), m->FractalDimension);
  201.        NewMountain = AllocMatrix(Float, NewSize+1, NewSize+1);
  202.        for (x = 0; x <= NewSize; x++)
  203.            for (y = 0; y <= NewSize; y++) {
  204.                XShort = (x+1)>>1;
  205.                XFar   = (x|1) - XShort;
  206.                YShort = (y+1)>>1;
  207.                YFar   = (y|1) - YShort;
  208.                NewMountain[x][y] =
  209.                    (OldMountain[XFar][YFar] +
  210.                     3 * (OldMountain[XFar][YShort] +
  211.                          OldMountain[XShort][YFar] +
  212.                          3 * OldMountain[XShort][YShort])) * 0.0625 +
  213.                     GAUSS(m->Number, x, y, i) * Scale;
  214.            }
  215.        FreeMatrix(OldMountain, OldSize+1);
  216.        OldMountain = NewMountain; OldSize = NewSize;
  217.    }
  218.    m->EvolvedData = OldMountain;
  219.    m->EvolvedSize = OldSize;
  220. }
  221.  
  222. /*
  223.  * Calculate how many cells on every iterationlevel
  224.  */
  225. static void MountainSetupLevelList(m)
  226. Mountain  *m;
  227. {
  228. int i, Ilevel, Elevel, Dlevel;
  229.  
  230.     Ilevel = m->IterationLevel;
  231.     Elevel = m->EvolveLevel;
  232.  
  233.     Dlevel = Ilevel - Elevel + 1;
  234.     m->MaxSize = (int *) share_malloc(Dlevel * sizeof(int));
  235.     m->MaxSize[0] = MyPower(2, Elevel) * (m->StartGridSize - 1) +  1;
  236.     for (i = 1; i < Dlevel; i++)
  237.          m->MaxSize[i] = 2 * m->MaxSize[i - 1] - 1;
  238. }
  239.  
  240. /*
  241.  * Allocate memory for cache on every iterationlevel
  242.  */
  243. static void MountainSetupCache(m)
  244. Mountain *m;
  245. {
  246. int   i,j, Ilevel, Elevel, Dlevel;
  247.  
  248.     Ilevel = m->IterationLevel;
  249.     Elevel = m->EvolveLevel;
  250.  
  251.     if (Dlevel = Ilevel - Elevel) {
  252.         m->Cache1 = (HeightCache **) share_malloc(Dlevel * sizeof(HeightCache *));
  253.         m->Cache2 = (HeightCache **) share_malloc(Dlevel * sizeof(HeightCache *));
  254.         for (i = 0; i < Dlevel; i++)  {
  255.              m->Cache1[i] = (HeightCache *) share_malloc((2 * m->MaxSize[i + 1] + 1) * sizeof(HeightCache));
  256.              m->Cache2[i] = (HeightCache *) share_malloc((2 * m->MaxSize[i + 1] + 1) * sizeof(HeightCache));
  257.              for (j = 0; j <= 2 * m->MaxSize[i + 1]; j++) {
  258.                   m->Cache1[i][j].x = -1;
  259.                   m->Cache2[i][j].x = -1;
  260.              }
  261.         }
  262.     }
  263. }
  264.  
  265.  
  266. /* 
  267.  * Setup list to keep track of which cells are needed
  268.  * on the different iterationlevels
  269.  */
  270. static void MountainSetupCalcList(m)
  271. Mountain   *m;
  272. {
  273. int       i, j;
  274. CalcList  *Chasingcl, *Cl;
  275.  
  276.     Chasingcl   = (CalcList *) share_malloc(sizeof(CalcList));
  277.     for (i = 0; i <= 2; i++)
  278.          for (j = 0; j <= 2; j++)
  279.               Chasingcl->ToDo[i][j] = FALSE;
  280.     m->ToCalcList       = Chasingcl;
  281.     Chasingcl->Previous = NULL;
  282.     for (i = m->IterationLevel; i > m->EvolveLevel + 1; i--) {
  283.         Cl              = (CalcList *) share_malloc(sizeof(CalcList));
  284.         Chasingcl->Next = Cl;
  285.         Cl->Previous    = Chasingcl;
  286.         Chasingcl       = Cl;
  287.     }
  288.     Chasingcl->Next = NULL;
  289. }
  290.  
  291.  
  292. /*
  293.  * Intersect ray with mountain.
  294.  */
  295. int MountainIntersect(mountain, ray, mindist, maxdist)
  296. Mountain   *mountain;
  297. Ray        *ray;
  298. Float       mindist, *maxdist;
  299. {
  300.         Vector  currpos;
  301.         Float   offset, currmaxdist;
  302.         int     x, y, sum;
  303.         int     hit=FALSE;
  304.         Float   Abs_t_X, Delta_t_X;
  305.         int     incX, boundX;
  306.         Float   Abs_t_Y, Delta_t_Y;
  307.         int     incY, boundY;
  308.         Float   Delta_w_X_z, Delta_w_Y_z;
  309.         Float   wx_z, wy_z, currpos_z;
  310.         int     triangleorder;
  311.         Float   Rx, Ry, Rsum;
  312.         
  313.  
  314.         MountainTests++;
  315.  
  316.         /*
  317.          * Find out where the projected ray hits the grid.
  318.          */
  319.         VecAddScaled(ray->pos, mindist, ray->dir, &currpos);
  320.         if (OutOf2DBounds(currpos, mountain->D2Bounds)) {
  321.              offset = *maxdist;
  322.              if (!MountainBoundsIntersect(ray, mountain->D2Bounds, mindist, &offset))
  323.                         return FALSE;
  324.              else {
  325.                   VecAddScaled(ray->pos, offset, ray->dir, &currpos);
  326.                   if (equal(currpos.x, mountain->D2Bounds[LOW][X]) ||
  327.                       equal(currpos.y, mountain->D2Bounds[LOW][Y]))
  328.                             triangleorder = 1;
  329.                   else triangleorder = -1;
  330.               }
  331.         } else {
  332.              offset = mindist;
  333.              Rx = x2grid(mountain, currpos.x);
  334.              Ry = y2grid(mountain, currpos.y);
  335.              Rsum = Rx + Ry;
  336.              sum = (int) Rx + (int) Ry;
  337.              if ((Rsum - sum) <= 0.5) triangleorder = 1;
  338.              else triangleorder = -1;
  339.         }
  340.  
  341.         x = x2grid(mountain, currpos.x);
  342.         if (x == mountain->GridMax) x--;
  343.  
  344.         if (fabs(ray->dir.x) < EPSILON) {
  345.                 Abs_t_X   = FAR_AWAY;
  346.                 Delta_t_X = 0;
  347.         } else if (ray->dir.x < 0.) {
  348.                 Abs_t_X     = offset + (grid2x(mountain, x) - currpos.x) / ray->dir.x;
  349.                 Delta_t_X   = mountain->XSize / (-ray->dir.x);
  350.                 incX = -1;  boundX = -1;
  351.         } else {
  352.                 Abs_t_X     = offset + (grid2x(mountain, x+1) - currpos.x) / ray->dir.x;
  353.                 Delta_t_X   = mountain->XSize / ray->dir.x;
  354.                 incX = 1;  boundX = mountain->GridMax;
  355.         }
  356.  
  357.         y = y2grid(mountain, currpos.y);
  358.         if (y == mountain->GridMax) y--;
  359.  
  360.         if (fabs(ray->dir.y) < EPSILON) {
  361.                 Abs_t_Y   = FAR_AWAY;
  362.                 Delta_t_Y = 0;
  363.          } else if (ray->dir.y < 0.) {
  364.                 Abs_t_Y     = offset + (grid2y(mountain, y) - currpos.y) / ray->dir.y;
  365.                 Delta_t_Y   = mountain->YSize / (-ray->dir.y);
  366.                 incY = -1;  boundY = -1;
  367.         } else {
  368.                 Abs_t_Y     = offset + (grid2y(mountain, y+1) - currpos.y) / ray->dir.y;
  369.                 Delta_t_Y   = mountain->YSize / ray->dir.y;
  370.                 incY = 1;  boundY = mountain->GridMax;
  371.         }
  372.  
  373.  
  374.         Delta_w_X_z = Delta_t_X * ray->dir.z;
  375.         Delta_w_Y_z = Delta_t_Y * ray->dir.z;
  376.  
  377.         wx_z = ray->pos.z + Abs_t_X * ray->dir.z;
  378.         wy_z = ray->pos.z + Abs_t_Y * ray->dir.z;
  379.  
  380.         currpos_z   = currpos.z;
  381.  
  382.         /*
  383.          * Traverse grid with DDA-algorithm
  384.          */
  385.         while(TRUE) {
  386.             if (Abs_t_X < Abs_t_Y) {
  387.                 currmaxdist = min(Abs_t_X, *maxdist);
  388.                 if (CheckGrid(mountain, x, y, ray, offset, &currmaxdist, triangleorder, currpos_z, wx_z)) {
  389.                     *maxdist = currmaxdist;
  390.                     return TRUE;
  391.                 }
  392.                 x += incX;
  393.                 if (x == boundX) break;
  394.                 triangleorder = incX;
  395.                 offset = currmaxdist;
  396.                 currpos_z  = wx_z;
  397.                 Abs_t_X   += Delta_t_X;
  398.                 wx_z      += Delta_w_X_z;
  399.             } else {
  400.                 currmaxdist = min(Abs_t_Y, *maxdist);
  401.                 if (CheckGrid(mountain, x, y, ray, offset, &currmaxdist, triangleorder, currpos_z, wy_z)) {
  402.                     *maxdist = currmaxdist;
  403.                     return TRUE;
  404.                 }
  405.                 y += incY;
  406.                 if (y == boundY) break;
  407.                 triangleorder = incY;
  408.                 offset = currmaxdist;
  409.                 currpos_z   = wy_z;
  410.                 Abs_t_Y   += Delta_t_Y;
  411.                 wy_z      += Delta_w_Y_z;
  412.            }
  413.         }
  414.         return FALSE;
  415. }
  416.  
  417.  
  418. int MountainBoundsIntersect(ray, D2Bounds, mindist, maxdist)
  419. Ray     *ray;
  420. Float    D2Bounds[2][2];
  421. Float    mindist;
  422. Float   *maxdist;
  423. {
  424.  
  425.     Float    t, tmin, tmax;
  426.     Float    dir, pos;
  427.  
  428.  
  429.     tmax = *maxdist;
  430.     tmin =  mindist;
  431.  
  432.     dir = ray->dir.x;
  433.     pos = ray->pos.x;
  434.  
  435.     if (dir < 0) {
  436.         t = (D2Bounds[LOW][X] - pos) / dir;
  437.         if (t < tmin) return FALSE;
  438.         if (t <= tmax) tmax = t;
  439.         t = (D2Bounds[HIGH][X] - pos) / dir;
  440.         if (t >= tmin) {
  441.             if (t > tmax) return FALSE;
  442.             tmin = t;
  443.         }
  444.     } else if (dir > 0) {
  445.         t = (D2Bounds[HIGH][X] - pos) / dir;
  446.         if (t < tmin) return FALSE;
  447.         if (t <= tmax) tmax = t;
  448.         t = (D2Bounds[LOW][X] - pos) / dir;
  449.         if (t >= tmin) {
  450.             if (t > tmax) return FALSE;
  451.             tmin = t;
  452.         }
  453.     } else if (pos < D2Bounds[LOW][X] ||
  454.                pos > D2Bounds[HIGH][X]) return FALSE;
  455.  
  456.     dir = ray->dir.y;
  457.     pos = ray->pos.y;
  458.  
  459.     if (dir < 0) {
  460.         t = (D2Bounds[LOW][Y] - pos) / dir;
  461.         if (t < tmin) return FALSE;
  462.         if (t <= tmax) tmax = t;
  463.         t = (D2Bounds[HIGH][Y] - pos) / dir;
  464.         if (t >= tmin) {
  465.             if (t > tmax) return FALSE;
  466.             tmin = t;
  467.         }
  468.     } else if (dir > 0) {
  469.         t = (D2Bounds[HIGH][Y] - pos) / dir;
  470.         if (t < tmin) return FALSE;
  471.         if (t <= tmax) tmax = t;
  472.         t = (D2Bounds[LOW][Y] - pos) / dir;
  473.         if (t >= tmin) {
  474.             if (t > tmax) return FALSE;
  475.             tmin = t;
  476.         }
  477.     } else if (pos < D2Bounds[LOW][Y] ||
  478.                pos > D2Bounds[HIGH][Y]) return FALSE;
  479.  
  480.  
  481.     if (tmin == mindist) {
  482.         if (tmax < *maxdist) {
  483.             *maxdist = tmax;
  484.             return TRUE;
  485.         }
  486.     } else {
  487.         if (tmin < *maxdist) {
  488.             *maxdist = tmin;
  489.             return TRUE;
  490.         }
  491.     }
  492.     return FALSE;
  493. }
  494.  
  495. #define CHECKCACHE(cache1, cache2, Rx, Ry, calclist, xx, yy)    \
  496.     for (ii = 0; ii <= 1; ii++, Rx++, Ry-=2) {                  \
  497.         for (jj = 0; jj <= 1; jj++, Ry++) {                     \
  498.             if ((cache1->x == Rx) && (cache1->y == Ry))         \
  499.                 calclist->Data[xx+ii][yy+jj] = cache1->z;       \
  500.             else if ((cache2->x == Rx) && (cache2->y == Ry))    \
  501.                 calclist->Data[xx+ii][yy+jj] = cache2->z;       \
  502.             else calclist->ToDo[xx+ii][yy+jj] = TRUE;           \
  503.             cache1++;                                           \
  504.             cache2--;                                           \
  505.         }                                                       \
  506.         cache1--;                                               \
  507.         cache2+=3;                                              \
  508.     }
  509.  
  510.  
  511. /*
  512.  * Calculate heightvalues at four corners of grid
  513.  * and check if ray possibly intersects the two triangles
  514.  * defined by the four points
  515.  */
  516. int CheckGrid(m, x, y, ray, offset, maxdist, triangleorder, z_near, z_far)
  517. Mountain    *m;
  518. int          x, y;
  519. Ray         *ray;
  520. Float        offset;
  521. Float       *maxdist;
  522. int          triangleorder;
  523. Float        z_near, z_far;
  524. {
  525.  
  526.    HeightCache *Cache1, *Cache2;
  527.    int          NX, NY, TX, TY, TTX, TTY, XP, YP;
  528.    int          level, i, j, xx, yy, ii, jj;
  529.    int          XShort, XFar, YShort, YFar;
  530.    CalcList    *TCL, *TCLP;
  531.    Vector       point1, point2, point3, edge1, edge2, edge3;
  532.    Float        scale;
  533.  
  534.    static CalcList  EvolveLevelData;
  535.  
  536.  
  537.  
  538.    level = m->IterationLevel - m->EvolveLevel - 1;
  539.    if (level >= 0) {
  540.        TCL = m->ToCalcList;
  541.        TCL->x = x;
  542.        TCL->y = y;
  543.        for (i = 0; i <= 1; i++)
  544.             for (j = 0; j<= 1; j++)
  545.                  TCL->ToDo[i][j] = FALSE;
  546.        /* 
  547.         * Check Cache if cell is available 
  548.         */
  549.        Cache1 = &m->Cache1[level][x+y];
  550.        Cache2 = &m->Cache2[level][m->MaxSize[level + 1] + x - y];
  551.        xx = x;  yy = y;
  552.        CHECKCACHE(Cache1, Cache2, xx, yy, TCL, 0, 0);
  553.        xx= x; yy=y;
  554.        level--;
  555.        /*  
  556.         * Calculate the coordinates of the cells necsarry for calculating
  557.         * the top cell, check cache if they are available 
  558.         */
  559.        for (TCLP = TCL, TCL = TCL->Next; TCL; TCLP = TCL, TCL = TCL->Next,
  560.             level--, xx = NX, yy = NY) {
  561.               NX = TCL->x = xx>>1;
  562.               NY = TCL->y = yy>>1;
  563.               for (i = 0; i <= 2; i++)
  564.                    for (j = 0; j <= 2; j++)
  565.                         TCL->ToDo[i][j] = FALSE;
  566.               for (i = 0; i <= 2; i++)
  567.                   for (j = 0; j <= 2; j++)
  568.                       if (TCLP->ToDo[i][j]) {
  569.                           TX = (xx+i)>>1;  TY = (yy+j)>>1;
  570.                           Cache1 = &m->Cache1[level][TX+TY];
  571.                           Cache2 = &m->Cache2[level][m->MaxSize[level + 1]+TX-TY];
  572.                           TTX = TX - NX;  TTY = TY - NY;
  573.                           CHECKCACHE(Cache1, Cache2, TX, TY, TCL, TTX, TTY);
  574.                       }
  575.        }
  576.        xx = xx>>1;  yy = yy>>1;
  577.        for (i = 0; i <= min(2, m->MaxSize[0] - xx); i++)
  578.            for (j = 0; j <= min(2, m->MaxSize[0] - yy); j++)
  579.                EvolveLevelData.Data[i][j] = m->EvolvedData[xx + i][yy + j];
  580.        level = 0;
  581.        XP    = xx; YP    = yy;
  582.        /*
  583.         * Calculate heightvalues of corners of cell who were
  584.         * not available in cache 
  585.         */
  586.        for (TCL = TCLP, TCLP = &EvolveLevelData; TCL; TCLP = TCL, TCL = TCL->Previous,
  587.             level++, XP = xx, YP = yy) {
  588.               xx = TCL->x;  yy = TCL->y;
  589.               Cache1 = &m->Cache1[level][xx + yy];
  590.               Cache2 = &m->Cache2[level][m->MaxSize[level + 1] + xx - yy];
  591.               scale = m->Scale * pow((m->StartGridSize/ (Float) m->MaxSize[level + 1]), m->FractalDimension);
  592.               for (i = 0; i<=2; i++) {
  593.                   for (j = 0; j <= 2; j++) {
  594.                       if (TCL->ToDo[i][j]) {
  595.                           XShort = (xx + 1 + i)>>1;
  596.                           XFar   = ((xx + i)|1) - XShort;
  597.                           YShort = (yy + 1 + j)>>1;
  598.                           YFar   = ((yy + j)|1) - YShort;
  599.                           XShort -= XP;
  600.                           XFar -= XP;
  601.                           YShort -= YP;
  602.                           YFar -= YP;
  603.                           TCL->Data[i][j] = (3*(TCLP->Data[XShort][YFar]
  604.                                                 + TCLP->Data[XFar][YShort]
  605.                                                 + 3 * TCLP->Data[XShort][YShort])
  606.                                              + TCLP->Data[XFar][YFar]) * 0.0625 +
  607.                                              GAUSS(m->Number, xx + i, yy + j, level + m->EvolveLevel + 1) * scale;
  608.  
  609.                           if (WhichCache[i][j]) {
  610.                               Cache1->x = xx + i;
  611.                               Cache1->y = yy + j;
  612.                               Cache1->z = TCL->Data[i][j];
  613.                           } else {
  614.                               Cache2->x = xx + i;
  615.                               Cache2->y = yy + j;
  616.                               Cache2->z = TCL->Data[i][j];
  617.                           }
  618.                       }
  619.                       Cache1++;
  620.                       Cache2--;
  621.                   }
  622.                   Cache1 -= 2;
  623.                   Cache2 += 4;
  624.               }
  625.           }
  626.    } else {
  627.        for (i = 0; i <= 1; i++)
  628.             for (j = 0; j <= 1; j++)
  629.                  m->ToCalcList->Data[i][j] = m->EvolvedData[x+i][y+j];
  630.    }
  631.    TCL = m->ToCalcList;
  632.    /*
  633.     * Check if ray crosses cell at right height 
  634.     */
  635.    if ((min(z_near, z_far) > max(max(max(TCL->Data[0][0], TCL->Data[0][1]), TCL->Data[1][0]), TCL->Data[1][1])) ||
  636.        (max(z_near, z_far) < min(min(min(TCL->Data[0][0], TCL->Data[0][1]), TCL->Data[1][0]), TCL->Data[1][1]))) return FALSE;
  637.    /*
  638.     * Test for ray-triangle intersection 
  639.     */
  640.    for (i = 1; i <=2 ; i++, triangleorder += 2)
  641.        if (triangleorder == 1) {
  642.            point1.x = (x * m->XSize) + m->D2Bounds[LOW][X];
  643.            point1.y = (y * m->YSize) + m->D2Bounds[LOW][Y];
  644.            point1.z = TCL->Data[0][0];
  645.            point2.x = ((x+1) * m->XSize) + m->D2Bounds[LOW][X];
  646.            point2.y = (y * m->YSize) + m->D2Bounds[LOW][Y];
  647.            point2.z = TCL->Data[1][0];
  648.            point3.x = (x * m->XSize) + m->D2Bounds[LOW][X];
  649.            point3.y = ((y+1) * m->YSize) + m->D2Bounds[LOW][Y];
  650.            point3.z = TCL->Data[0][1];
  651.            if (IntersectTriangle(m, ray, &point1, &point2, &point3, offset, maxdist)) return TRUE;
  652.        } else {
  653.            point1.x = ((x+1) * m->XSize) + m->D2Bounds[LOW][X];
  654.            point1.y = ((y+1) * m->YSize) + m->D2Bounds[LOW][Y];
  655.            point1.z = TCL->Data[1][1];
  656.            point2.x = (x * m->XSize) + m->D2Bounds[LOW][X];
  657.            point2.y = ((y+1) * m->YSize) + m->D2Bounds[LOW][Y];
  658.            point2.z = TCL->Data[0][1];
  659.            point3.x = ((x+1) * m->XSize) + m->D2Bounds[LOW][X];
  660.            point3.y = (y * m->YSize) + m->D2Bounds[LOW][Y];
  661.            point3.z = TCL->Data[1][0];
  662.            if (IntersectTriangle(m, ray, &point1, &point2, &point3, offset, maxdist)) return TRUE;
  663.        }
  664.    return FALSE;
  665. }
  666.  
  667. /* 
  668.  * Check if ray intersects triangle 
  669.  */
  670.  
  671. int IntersectTriangle(m, ray, point1, point2, point3, mindist, maxdist)
  672. Mountain   *m;
  673. Ray        *ray;
  674. Vector     *point1, *point2, *point3;
  675. Float       mindist, *maxdist;
  676. {
  677.     Vector  edge1, edge2, edge3;
  678.     Float   k, s;
  679.     Vector  normal, absnorm;
  680.     Vector  pos, dir;
  681.     Float   qi1, qi2, b0, b1, b2;
  682.  
  683.  
  684.     VecSub(*point2, *point1, &edge1);
  685.     VecSub(*point3, *point2, &edge2);
  686.     VecSub(*point1, *point3, &edge3);
  687.  
  688.     /* 
  689.      * Find plane normal. 
  690.      */
  691.     VecCross(&edge1, &edge2, &normal);
  692.  
  693.     pos = ray->pos;
  694.     dir = ray->dir;
  695.     /*
  696.      * Plane intersection.
  697.      */
  698.     k = dotp(&normal, &dir);
  699.     if (fabs(k) < EPSILON)
  700.         return FALSE;
  701.     s = (dotp(&normal, point1) - dotp(&normal, &pos)) / k;
  702.     if (s < mindist || s > *maxdist)
  703.                 return FALSE;
  704.  
  705.     /*
  706.      * Find "dominant" part of normal vector.
  707.      */
  708.     absnorm.x = fabs(normal.x);
  709.     absnorm.y = fabs(normal.y);
  710.     absnorm.z = fabs(normal.z);
  711.  
  712.     if (absnorm.x > absnorm.y && absnorm.x > absnorm.z) {
  713.         qi1 = pos.y + s * dir.y;
  714.         qi2 = pos.z + s * dir.z;
  715.         b0 = (edge2.y * (qi2 - point2->z) -
  716.               edge2.z * (qi1 - point2->y)) / normal.x;
  717.         if (b0 < 0. || b0 > 1.) return FALSE;
  718.         b1 = (edge3.y * (qi2 - point3->z) -
  719.               edge3.z * (qi1 - point3->y)) / normal.x;
  720.         if (b1 < 0. || b1 > 1.) return FALSE;
  721.         b2 = (edge1.y * (qi2 - point1->z) -
  722.               edge1.z * (qi1 - point1->y)) / normal.x;
  723.         if (b2 < 0. || b2 > 1.) return FALSE;
  724.     } else if (absnorm.y > absnorm.z) {
  725.         qi1 = pos.x + s * dir.x;
  726.         qi2 = pos.z + s * dir.z;
  727.         b0 = (edge2.z * (qi1 - point2->x) -
  728.               edge2.x * (qi2 - point2->z)) / normal.y;
  729.         if (b0 < 0. || b0 > 1.) return FALSE;
  730.         b1 = (edge3.z * (qi1 - point3->x) -
  731.               edge3.x * (qi2 - point3->z)) / normal.y;
  732.         if (b1 < 0. || b1 > 1.) return FALSE;
  733.         b2 = (edge1.z * (qi1 - point1->x) -
  734.               edge1.x * (qi2 - point1->z)) / normal.y;
  735.         if (b2 < 0. || b2 > 1.) return FALSE;
  736.     } else {
  737.         qi1 = pos.x + s * dir.x;
  738.         qi2 = pos.y + s * dir.y;
  739.         b0 = (edge2.x * (qi2 - point2->y) -
  740.               edge2.y * (qi1 - point2->x)) / normal.z;
  741.         if (b0 < 0. || b0 > 1.) return FALSE;
  742.         b1 = (edge3.x * (qi2 - point3->y) -
  743.               edge3.y * (qi1 - point3->x)) / normal.z;
  744.         if (b1 < 0. || b1 > 1.) return FALSE;
  745.         b2 = (edge1.x * (qi2 - point1->y) -
  746.               edge1.y * (qi1 - point1->x)) / normal.z;
  747.         if (b2 < 0. || b2 > 1.) return FALSE;
  748.     }
  749.     MountainHits++;
  750.     *maxdist = s;
  751.     m->Normal = normal;
  752.     VecAddScaled(ray->pos, s, ray->dir, &m->IntersectPos);
  753.  
  754.     return TRUE;
  755. }
  756.  
  757. #define VecEqual(a,b) ((equal((a).x, (b).x)) && (equal((a).y, (b).y)) && (equal((a).z, (b).z)))
  758.  
  759. /*
  760.  * Calculate normal on triangle in pos
  761.  */
  762. int MountainNormal(m, pos, nrm, gnrm)
  763. Mountain *m;
  764. Vector   *pos, *nrm, *gnrm;
  765. {
  766.     if (!VecEqual(*pos, m->IntersectPos))
  767.         RLerror(RL_ABORT, "Can't compute mountain normal in point other than the last intersection.\n");
  768.         *gnrm = m->Normal;
  769.  
  770.         VecNormalize(gnrm);
  771.         *nrm = *gnrm;
  772.         return FALSE;
  773. }
  774.  
  775.  
  776. char *MountainName()
  777. {
  778.         return mountainName;
  779. }
  780.  
  781. void MountainStats(tests, hits)
  782. unsigned long *tests, *hits;
  783. {
  784.         *tests = MountainTests;
  785.         *hits = MountainHits;
  786. }
  787.  
  788. void MountainUV(m, pos, norm, uv, dpdu, dpdv)
  789. Mountain   *m;
  790. Vector     *pos, *norm, *dpdu, *dpdv;
  791. Vec2d      *uv;
  792. {
  793.         uv->u = pos->x;
  794.         uv->v = pos->y;
  795.         if (dpdu) {
  796.               dpdu->x = 1.;
  797.               dpdu->y = dpdu->z = 0.;
  798.               dpdv->x = dpdv->z = 0.;
  799.               dpdv->y = 1.;
  800.         }
  801. }
  802.  
  803. void MountainBounds(m, bounds)
  804. Mountain *m;
  805. Float     bounds[2][3];
  806. {
  807.         bounds[LOW][X] = m->D2Bounds[LOW][X];
  808.         bounds[LOW][Y] = m->D2Bounds[LOW][Y];
  809.         bounds[HIGH][X] = m->D2Bounds[HIGH][X];
  810.         bounds[HIGH][Y] = m->D2Bounds[HIGH][Y];
  811.         bounds[LOW][Z] = -FAR_AWAY;
  812.         bounds[HIGH][Z] = FAR_AWAY;
  813. }
  814.  
  815. void MountainMethodRegister(meth)
  816. UserMethodType  meth;
  817. {
  818.         if (iMountainMethods)
  819.             iMountainMethods->user = meth;
  820. }
  821.  
  822.  
  823. Methods *MountainMethods()
  824. {
  825.         if (iMountainMethods == (Methods *)NULL) {
  826.                 iMountainMethods = MethodsCreate();
  827.                 iMountainMethods->create = (GeomCreateFunc *)MountainCreate;
  828.                 iMountainMethods->methods = MountainMethods;
  829.                 iMountainMethods->name = MountainName;
  830.                 iMountainMethods->intersect = MountainIntersect;
  831.                 iMountainMethods->normal = MountainNormal;
  832.                 iMountainMethods->uv = MountainUV;
  833.                 iMountainMethods->bounds = MountainBounds;
  834.                 iMountainMethods->stats = MountainStats;
  835.                 iMountainMethods->checkbounds = FALSE;
  836.                 iMountainMethods->closed = FALSE;
  837.         }
  838.         return iMountainMethods;
  839. }
  840.